#############################################
## PROGRAM NAME: 204a_checks_msa.R         ##
## AUTHOR: MATT MLECZKO                    ##
## INPUTS:                                 ##
##         202a_m1_msa.Rda -               ##
##         202a_m12_msa.Rda                ##
##                                         ##
## OUTPUTS:                                ##
##         204a_wts_out_msa.csv            ##
##                                         ## 
## PURPOSE: MSA MSM model checks           ##
##                                         ##
## LIST OF UPDATES:                        ##
#############################################

#log <- file(# USER DEFINED PATH AND FILE NAME HERE #)
#sink(log, append=TRUE)
#sink(log, append=TRUE, type="message")

## load libraries ##
library(tidyverse)

data_path <- # USER DEFINED PATH HERE #
setwd(data_path)

## load the results ##

load("202a_m1_msa.Rda")
load("202a_m2_msa.Rda")
load("202a_m3_msa.Rda")
load("202a_m4_msa.Rda")
load("202a_m5_msa.Rda")
load("202a_m6_msa.Rda")
load("202a_m7_msa.Rda")
load("202a_m8_msa.Rda")
load("202a_m9_msa.Rda")
load("202a_m10_msa.Rda")
load("202a_m11_msa.Rda")
load("202a_m12_msa.Rda")


#############
## weights ##
#############

ex.wts <- function(i, inlist) {
  
  res <- inlist[[i]][2]
  res.new <- as.vector(res)
  
  return(res.new)
  
}

w1.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res1)

w1 <- unlist(w1.res)
summary(w1)

w2.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res2)

w2 <- unlist(w2.res)
summary(w2)

w3.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res3)

w3 <- unlist(w3.res)
summary(w3)

w4.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res4)

w4 <- unlist(w4.res)
summary(w4)

w5.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res5)

w5 <- unlist(w5.res)
summary(w5)

w6.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res6)

w6 <- unlist(w6.res)
summary(w6)

w7.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res7)

w7 <- unlist(w7.res)
summary(w7)

w8.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res8)

w8 <- unlist(w8.res)
summary(w8)

w9.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res9)

w9 <- unlist(w9.res)
summary(w9)

w10.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res10)

w10 <- unlist(w10.res)
summary(w10)


w11.res <- sapply(seq(1,100,by=1),
                 ex.wts,
                 inlist = msa.res11)

w11 <- unlist(w11.res)
summary(w11)

w12.res <- sapply(seq(1,100,by=1),
                  ex.wts,
                  inlist = msa.res12)

w12 <- unlist(w12.res)
summary(w12)


## put all wts together ##

wts.out <- data.frame(min = c(summary(w1)[1],
                              summary(w2)[1],
                              summary(w3)[1],
                              summary(w4)[1],
                              summary(w5)[1],
                              summary(w6)[1],
                              summary(w7)[1],
                              summary(w8)[1],
                              summary(w9)[1],
                              summary(w10)[1],
                              summary(w11)[1],
                              summary(w12)[1]),
                      q1 = c(summary(w1)[2],
                             summary(w2)[2],
                             summary(w3)[2],
                             summary(w4)[2],
                             summary(w5)[2],
                             summary(w6)[2],
                             summary(w7)[2],
                             summary(w8)[2],
                             summary(w9)[2],
                             summary(w10)[2],
                             summary(w11)[2],
                             summary(w12)[2]),
                      median = c(summary(w1)[3],
                                 summary(w2)[3],
                                 summary(w3)[3],
                                 summary(w4)[3],
                                 summary(w5)[3],
                                 summary(w6)[3],
                                 summary(w7)[3],
                                 summary(w8)[3],
                                 summary(w9)[3],
                                 summary(w10)[3],
                                 summary(w11)[3],
                                 summary(w12)[3]),
                      mean = c(summary(w1)[4],
                               summary(w2)[4],
                               summary(w3)[4],
                               summary(w4)[4],
                               summary(w5)[4],
                               summary(w6)[4],
                               summary(w7)[4],
                               summary(w8)[4],
                               summary(w9)[4],
                               summary(w10)[4],
                               summary(w11)[4],
                               summary(w12)[4]),
                      q3 = c(summary(w1)[5],
                             summary(w2)[5],
                             summary(w3)[5],
                             summary(w4)[5],
                             summary(w5)[5],
                             summary(w6)[5],
                             summary(w7)[5],
                             summary(w8)[5],
                             summary(w9)[5],
                             summary(w10)[5],
                             summary(w11)[5],
                             summary(w12)[5]),
                      max = c(summary(w1)[6],
                              summary(w2)[6],
                              summary(w3)[6],
                              summary(w4)[6],
                              summary(w5)[6],
                              summary(w6)[6],
                              summary(w7)[6],
                              summary(w8)[6],
                              summary(w9)[6],
                              summary(w10)[6],
                              summary(w11)[6],
                              summary(w12)[6]),
                      sd = c(sd(w1,na.rm=T),
                             sd(w2,na.rm=T),
                             sd(w3,na.rm=T),
                             sd(w4,na.rm=T),
                             sd(w5,na.rm=T),
                             sd(w6,na.rm=T),
                             sd(w7,na.rm=T),
                             sd(w8,na.rm=T),
                             sd(w9,na.rm=T),
                             sd(w10,na.rm=T),
                             sd(w11,na.rm=T),
                             sd(w12,na.rm=T)))

write.csv(wts.out, 
          "204a_wts_out_msa.csv")

####################
## balance checks ##
####################

unwt <- function(i, k, inlist) {
  
  out.bal <- list()
  vars.unwt.st <- c()
  vars.wt.st <- c()
  
  for(j in 1:15){
    
    ## outermost ##
    t2 <- inlist[[i]][k]
    ## var j ##
    t2v <- t2[[1]][j]
    #print(t2v)
    ## unweighted ##
    t2v.uw <- t2v[[1]][1]
    #print(t2v.uw)
    ## unwt est ## 
    vars.unwt.st[j] <- t2v.uw[[1]]$estimate[2]/t2v.uw[[1]]$std.error[2]
    ## weighted ##
    t2v.w <- t2v[[1]][2]
    ## wt est ## 
    vars.wt.st[j] <- t2v.w[[1]]$estimate[2]/t2v.w[[1]]$std.error[2]
    
  }
  
  return(vars.unwt.st)

}


wt <- function(i, k, inlist) {
  
  out.bal <- list()
  vars.unwt.st <- c()
  vars.wt.st <- c()
  
  for(j in 1:15){
    
    ## outermost ##
    t2 <- inlist[[i]][k]
    ## var j ##
    t2v <- t2[[1]][j]
    ## unweighted ##
    t2v.uw <- t2v[[1]][1]
    ## unwt est ## 
    vars.unwt.st[j] <- t2v.uw[[1]]$estimate[2]/t2v.uw[[1]]$std.error[2]
    ## weighted ##
    t2v.w <- t2v[[1]][2]
    ## wt est ## 
    vars.wt.st[j] <- t2v.w[[1]]$estimate[2]/t2v.w[[1]]$std.error[2]
    
  }
  
  
  return(vars.wt.st)
  
  
}

## M1 - RES ##

m1.t2.unwt.res <- sapply(seq(1,100,by=1),
                      unwt,
                      k=6,
                      inlist = msa.res1)
  

m1.t2.unwt.avg <- rowMeans(bind_cols(m1.t2.unwt.res))
  
m1.t2.wt.res <- sapply(seq(1,100,by=1),
                    wt,
                    k=6,
                    inlist = msa.res1)
    
m1.t2.wt.avg <- rowMeans(bind_cols(m1.t2.wt.res))
  
m1.t2.m1.bal <- as.data.frame(cbind(m1.t2.unwt.avg,
                                    m1.t2.wt.avg)) %>%
    mutate(var = c("pop_total","pop_density",
                   "per_18t34","per_35t64",
                   "per_65over","entropy_er",
                   "per_child", "per_nomove",
                   "per_fb","hownrt",
                   "median_yr_str","median_hhld_inc",
                   "hhld_pov_rt","unemr","gini")) %>%
    select(var,
           m1.t2.unwt.avg,
           m1.t2.wt.avg)
  
  
m1.t1.unwt.res <- sapply(seq(1,100,by=1),
                      unwt,
                      k=7,
                      inlist = msa.res1)
  
m1.t1.unwt.avg <- rowMeans(bind_cols(m1.t1.unwt.res))
  
m1.t1.wt.res <- sapply(seq(1,100,by=1),
                    wt,
                    k=7,
                    inlist = msa.res1)
  
m1.t1.wt.avg <- rowMeans(bind_cols(m1.t1.wt.res))
  
m1.t1.m1.bal <- as.data.frame(cbind(m1.t1.unwt.avg,
                                   m1.t1.wt.avg)) %>%
    mutate(var = c("pop_total","pop_density",
                   "per_18t34","per_35t64",
                   "per_65over","entropy_er",
                   "per_child", "per_nomove",
                   "per_fb","hownrt",
                   "median_yr_str","median_hhld_inc",
                   "hhld_pov_rt","unemr","gini")) %>%
    select(var,
           m1.t1.unwt.avg,
           m1.t1.wt.avg)
  
m1.t1t2.m1.bal <- m1.t2.m1.bal %>%
    inner_join(m1.t1.m1.bal, "var")
  


## M2 - RES ##

m2.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res2)


m2.t2.unwt.avg <- rowMeans(bind_cols(m2.t2.unwt.res))

m2.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res2)

m2.t2.wt.avg <- rowMeans(bind_cols(m2.t2.wt.res))

m2.t2.m2.bal <- as.data.frame(cbind(m2.t2.unwt.avg,
                                    m2.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m2.t2.unwt.avg,
         m2.t2.wt.avg)


m2.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res2)

m2.t1.unwt.avg <- rowMeans(bind_cols(m2.t1.unwt.res))

m2.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res2)

m2.t1.wt.avg <- rowMeans(bind_cols(m2.t1.wt.res))

m2.t1.m2.bal <- as.data.frame(cbind(m2.t1.unwt.avg,
                                    m2.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m2.t1.unwt.avg,
         m2.t1.wt.avg)

m2.t1t2.m2.bal <- m2.t2.m2.bal %>%
  inner_join(m2.t1.m2.bal, "var")



## M3 - RES ##

m3.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res3)


m3.t2.unwt.avg <- rowMeans(bind_cols(m3.t2.unwt.res))

m3.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res3)

m3.t2.wt.avg <- rowMeans(bind_cols(m3.t2.wt.res))

m3.t2.m3.bal <- as.data.frame(cbind(m3.t2.unwt.avg,
                                    m3.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m3.t2.unwt.avg,
         m3.t2.wt.avg)


m3.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res3)

m3.t1.unwt.avg <- rowMeans(bind_cols(m3.t1.unwt.res))

m3.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res3)

m3.t1.wt.avg <- rowMeans(bind_cols(m3.t1.wt.res))

m3.t1.m3.bal <- as.data.frame(cbind(m3.t1.unwt.avg,
                                    m3.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m3.t1.unwt.avg,
         m3.t1.wt.avg)

m3.t1t2.m3.bal <- m3.t2.m3.bal %>%
  inner_join(m3.t1.m3.bal, "var")


## M4 - RES ##

m4.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res4)


m4.t2.unwt.avg <- rowMeans(bind_cols(m4.t2.unwt.res))

m4.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res4)

m4.t2.wt.avg <- rowMeans(bind_cols(m4.t2.wt.res))

m4.t2.m4.bal <- as.data.frame(cbind(m4.t2.unwt.avg,
                                    m4.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m4.t2.unwt.avg,
         m4.t2.wt.avg)


m4.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res4)

m4.t1.unwt.avg <- rowMeans(bind_cols(m4.t1.unwt.res))

m4.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res4)

m4.t1.wt.avg <- rowMeans(bind_cols(m4.t1.wt.res))

m4.t1.m4.bal <- as.data.frame(cbind(m4.t1.unwt.avg,
                                    m4.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m4.t1.unwt.avg,
         m4.t1.wt.avg)

m4.t1t2.m4.bal <- m4.t2.m4.bal %>%
  inner_join(m4.t1.m4.bal, "var")




## M5 - RES ##

m5.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res5)


m5.t2.unwt.avg <- rowMeans(bind_cols(m5.t2.unwt.res))

m5.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res5)

m5.t2.wt.avg <- rowMeans(bind_cols(m5.t2.wt.res))

m5.t2.m5.bal <- as.data.frame(cbind(m5.t2.unwt.avg,
                                    m5.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m5.t2.unwt.avg,
         m5.t2.wt.avg)


m5.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res5)

m5.t1.unwt.avg <- rowMeans(bind_cols(m5.t1.unwt.res))

m5.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res5)

m5.t1.wt.avg <- rowMeans(bind_cols(m5.t1.wt.res))

m5.t1.m5.bal <- as.data.frame(cbind(m5.t1.unwt.avg,
                                    m5.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m5.t1.unwt.avg,
         m5.t1.wt.avg)

m5.t1t2.m5.bal <- m5.t2.m5.bal %>%
  inner_join(m5.t1.m5.bal, "var")




## M6 - RES ##

m6.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res6)


m6.t2.unwt.avg <- rowMeans(bind_cols(m6.t2.unwt.res))

m6.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res6)

m6.t2.wt.avg <- rowMeans(bind_cols(m6.t2.wt.res))

m6.t2.m6.bal <- as.data.frame(cbind(m6.t2.unwt.avg,
                                    m6.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m6.t2.unwt.avg,
         m6.t2.wt.avg)


m6.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res6)

m6.t1.unwt.avg <- rowMeans(bind_cols(m6.t1.unwt.res))

m6.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res6)

m6.t1.wt.avg <- rowMeans(bind_cols(m6.t1.wt.res))

m6.t1.m6.bal <- as.data.frame(cbind(m6.t1.unwt.avg,
                                    m6.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m6.t1.unwt.avg,
         m6.t1.wt.avg)

m6.t1t2.m6.bal <- m6.t2.m6.bal %>%
  inner_join(m6.t1.m6.bal, "var")



## M7 - RES ##

m7.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res7)


m7.t2.unwt.avg <- rowMeans(bind_cols(m7.t2.unwt.res))

m7.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res7)

m7.t2.wt.avg <- rowMeans(bind_cols(m7.t2.wt.res))

m7.t2.m7.bal <- as.data.frame(cbind(m7.t2.unwt.avg,
                                    m7.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m7.t2.unwt.avg,
         m7.t2.wt.avg)


m7.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res7)

m7.t1.unwt.avg <- rowMeans(bind_cols(m7.t1.unwt.res))

m7.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res7)

m7.t1.wt.avg <- rowMeans(bind_cols(m7.t1.wt.res))

m7.t1.m7.bal <- as.data.frame(cbind(m7.t1.unwt.avg,
                                    m7.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m7.t1.unwt.avg,
         m7.t1.wt.avg)

m7.t1t2.m7.bal <- m7.t2.m7.bal %>%
  inner_join(m7.t1.m7.bal, "var")


## M8 - RES ##

m8.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res8)


m8.t2.unwt.avg <- rowMeans(bind_cols(m8.t2.unwt.res))

m8.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res8)

m8.t2.wt.avg <- rowMeans(bind_cols(m8.t2.wt.res))

m8.t2.m8.bal <- as.data.frame(cbind(m8.t2.unwt.avg,
                                    m8.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m8.t2.unwt.avg,
         m8.t2.wt.avg)


m8.t1.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=7,
                         inlist = msa.res8)

m8.t1.unwt.avg <- rowMeans(bind_cols(m8.t1.unwt.res))

m8.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res8)

m8.t1.wt.avg <- rowMeans(bind_cols(m8.t1.wt.res))

m8.t1.m8.bal <- as.data.frame(cbind(m8.t1.unwt.avg,
                                    m8.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m8.t1.unwt.avg,
         m8.t1.wt.avg)

m8.t1t2.m8.bal <- m8.t2.m8.bal %>%
  inner_join(m8.t1.m8.bal, "var")


## M9 - RES ##

m9.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res9)


m9.t2.unwt.avg <- rowMeans(bind_cols(m9.t2.unwt.res))

m9.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res9)

m9.t2.wt.avg <- rowMeans(bind_cols(m9.t2.wt.res))

m9.t2.m9.bal <- as.data.frame(cbind(m9.t2.unwt.avg,
                                    m9.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m9.t2.unwt.avg,
         m9.t2.wt.avg)


m9.t1.unwt.res <- sapply(seq(1,100,by=9),
                         unwt,
                         k=7,
                         inlist = msa.res9)

m9.t1.unwt.avg <- rowMeans(bind_cols(m9.t1.unwt.res))

m9.t1.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=7,
                       inlist = msa.res9)

m9.t1.wt.avg <- rowMeans(bind_cols(m9.t1.wt.res))

m9.t1.m9.bal <- as.data.frame(cbind(m9.t1.unwt.avg,
                                    m9.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m9.t1.unwt.avg,
         m9.t1.wt.avg)

m9.t1t2.m9.bal <- m9.t2.m9.bal %>%
  inner_join(m9.t1.m9.bal, "var")



## M10 - RES ##

m10.t2.unwt.res <- sapply(seq(1,100,by=1),
                         unwt,
                         k=6,
                         inlist = msa.res10)


m10.t2.unwt.avg <- rowMeans(bind_cols(m10.t2.unwt.res))

m10.t2.wt.res <- sapply(seq(1,100,by=1),
                       wt,
                       k=6,
                       inlist = msa.res10)

m10.t2.wt.avg <- rowMeans(bind_cols(m10.t2.wt.res))

m10.t2.m10.bal <- as.data.frame(cbind(m10.t2.unwt.avg,
                                    m10.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m10.t2.unwt.avg,
         m10.t2.wt.avg)


m10.t1.unwt.res <- sapply(seq(1,100,by=10),
                         unwt,
                         k=7,
                         inlist = msa.res10)

m10.t1.unwt.avg <- rowMeans(bind_cols(m10.t1.unwt.res))

m10.t1.wt.res <- sapply(seq(1,100,by=10),
                       wt,
                       k=7,
                       inlist = msa.res10)

m10.t1.wt.avg <- rowMeans(bind_cols(m10.t1.wt.res))

m10.t1.m10.bal <- as.data.frame(cbind(m10.t1.unwt.avg,
                                    m10.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m10.t1.unwt.avg,
         m10.t1.wt.avg)

m10.t1t2.m10.bal <- m10.t2.m10.bal %>%
  inner_join(m10.t1.m10.bal, "var")



## M11 - RES ##

m11.t2.unwt.res <- sapply(seq(1,100,by=1),
                          unwt,
                          k=6,
                          inlist = msa.res11)


m11.t2.unwt.avg <- rowMeans(bind_cols(m11.t2.unwt.res))

m11.t2.wt.res <- sapply(seq(1,100,by=1),
                        wt,
                        k=6,
                        inlist = msa.res11)

m11.t2.wt.avg <- rowMeans(bind_cols(m11.t2.wt.res))

m11.t2.m11.bal <- as.data.frame(cbind(m11.t2.unwt.avg,
                                      m11.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m11.t2.unwt.avg,
         m11.t2.wt.avg)


m11.t1.unwt.res <- sapply(seq(1,100,by=1),
                          unwt,
                          k=7,
                          inlist = msa.res11)

m11.t1.unwt.avg <- rowMeans(bind_cols(m11.t1.unwt.res))

m11.t1.wt.res <- sapply(seq(1,100,by=1),
                        wt,
                        k=7,
                        inlist = msa.res11)

m11.t1.wt.avg <- rowMeans(bind_cols(m11.t1.wt.res))

m11.t1.m11.bal <- as.data.frame(cbind(m11.t1.unwt.avg,
                                      m11.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m11.t1.unwt.avg,
         m11.t1.wt.avg)

m11.t1t2.m11.bal <- m11.t2.m11.bal %>%
  inner_join(m11.t1.m11.bal, "var")



## M12 - RES ##

m12.t2.unwt.res <- sapply(seq(1,100,by=1),
                          unwt,
                          k=6,
                          inlist = msa.res12)


m12.t2.unwt.avg <- rowMeans(bind_cols(m12.t2.unwt.res))

m12.t2.wt.res <- sapply(seq(1,100,by=1),
                        wt,
                        k=6,
                        inlist = msa.res12)

m12.t2.wt.avg <- rowMeans(bind_cols(m12.t2.wt.res))

m12.t2.m12.bal <- as.data.frame(cbind(m12.t2.unwt.avg,
                                      m12.t2.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m12.t2.unwt.avg,
         m12.t2.wt.avg)


m12.t1.unwt.res <- sapply(seq(1,100,by=1),
                          unwt,
                          k=7,
                          inlist = msa.res12)

m12.t1.unwt.avg <- rowMeans(bind_cols(m12.t1.unwt.res))

m12.t1.wt.res <- sapply(seq(1,100,by=1),
                        wt,
                        k=7,
                        inlist = msa.res12)

m12.t1.wt.avg <- rowMeans(bind_cols(m12.t1.wt.res))

m12.t1.m12.bal <- as.data.frame(cbind(m12.t1.unwt.avg,
                                      m12.t1.wt.avg)) %>%
  mutate(var = c("pop_total","pop_density",
                 "per_18t34","per_35t64",
                 "per_65over","entropy_er",
                 "per_child", "per_nomove",
                 "per_fb","hownrt",
                 "median_yr_str","median_hhld_inc",
                 "hhld_pov_rt","unemr","gini")) %>%
  select(var,
         m12.t1.unwt.avg,
         m12.t1.wt.avg)

m12.t1t2.m12.bal <- m12.t2.m12.bal %>%
  inner_join(m12.t1.m12.bal, "var")


names(m1.t1t2.m1.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m2.t1t2.m2.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m3.t1t2.m3.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m4.t1t2.m4.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m5.t1t2.m5.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m6.t1t2.m6.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m7.t1t2.m7.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m8.t1t2.m8.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m9.t1t2.m9.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m10.t1t2.m10.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m11.t1t2.m11.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

names(m12.t1t2.m12.bal) <- c("var",
                           "t2.unwt.avg",
                           "t2.wt.avg",
                           "t1.unwt.avg",
                           "t1.wt.avg")

bal.rd <- rbind(m1.t1t2.m1.bal,
                m2.t1t2.m2.bal,
                m3.t1t2.m3.bal,
                m4.t1t2.m4.bal,
                m5.t1t2.m5.bal,
                m6.t1t2.m6.bal,
                m7.t1t2.m7.bal,
                m8.t1t2.m8.bal,
                m9.t1t2.m9.bal,
                m10.t1t2.m10.bal,
                m11.t1t2.m11.bal,
                m12.t1t2.m12.bal)

## plot ## 

par(mfrow=c(2,2))

plot(bal.rd$t2.wt.avg, col = "blue",
     main = "T2 Weighted",
     ylim = c(-5,5),
     ylab = "",
     xlab = "",
     xaxt='n')
abline(h = c(-2,2), col = "red")

plot(bal.rd$t2.unwt.avg, col = "blue",
     main = "T2 Unweighted",
     ylim = c(-5,5),
     ylab = "",
     xlab = "",
     xaxt='n')
abline(h = c(-2,2), col = "red")

plot(bal.rd$t1.wt.avg, col = "blue",
     main = "T1 weighted",
     ylim = c(-5,5),
     ylab = "",
     xlab = "",
     xaxt='n')
abline(h = c(-2,2), col = "red")

plot(bal.rd$t1.unwt.avg, col = "blue",
     main = "T1 Unweighted",
     ylim = c(-5,5),
     ylab = "",
     xlab = "",
     xaxt='n')
abline(h = c(-2,2), col = "red")


sum(bal.rd$t2.wt.avg <= 2 &
      bal.rd$t2.wt.avg >= -2)/nrow(bal.rd)

sum(bal.rd$t2.unwt.avg <= 2 &
      bal.rd$t2.unwt.avg >= -2)/nrow(bal.rd)

sum(bal.rd$t1.wt.avg <= 2 &
      bal.rd$t1.wt.avg >= -2)/nrow(bal.rd)

sum(bal.rd$t1.unwt.avg <= 2 &
      bal.rd$t1.unwt.avg >= -2)/nrow(bal.rd)

## END OF PROGRAM ##

#sink() 

